Fit density models to cod and flounder and check spatial overlap
Max Lindmark
05 December, 2022
Intro
In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.
Load packages & source functions
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(janitor)
library(forcats)
library(gapminder)
library(viridis)
library(ggridges)
library(raster)
library(icesDatras)
library(ggalluvial)
library(ggrepel)
library(ncdf4)
library(chron)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
library(quantreg)
#> Warning in .recacheSubclasses(def@className, def, env): undefined subclass
#> "numericVector" of class "Mnumeric"; definition not updated
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
library(sdmTMB)
library(sp)
library(raster)
library(ggcorrplot)
library(ggh4x)
library(visreg)
options(mc.cores = parallel::detectCores())
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")
theme_set(theme_plot())
# Continuous colors
options(ggplot2.continuous.colour = "viridis")
# Discrete colors
scale_colour_discrete <- function(...) {
scale_colour_brewer(palette = "Dark2")
}
scale_fill_discrete <- function(...) {
scale_fill_brewer(palette = "Dark2")
}Read and summarise data
cpue_full <- readr::read_csv("data/clean/catch_by_length_q1_q4.csv") %>%
janitor::clean_names() %>%
rename(X = x,
Y = y)
#> Rows: 1759609 Columns: 22
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): Country, haul.id, IDx, ices_rect, id_haul_stomach, species, haul.i...
#> dbl (14): density, year, lat, lon, quarter, Month, sub_div, length_cm, depth...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
# Summaries density by species group an haul
cpue <- cpue_full %>%
mutate(species2 = species,
species2 = ifelse(species == "cod" & length_cm >= 25, "large_cod", species),
species2 = ifelse(species == "cod" & length_cm < 25, "small_cod", species2)) %>%
group_by(haul_id, species2) %>%
summarise(density = sum(density)) %>%
ungroup()
#> mutate: new variable 'species2' (character) with 3 unique values and 0% NA
#> group_by: 2 grouping variables (haul_id, species2)
#> summarise: now 27,897 rows and 3 columns, one group variable remaining (haul_id)
#> ungroup: no grouping variables
# Trim data with quantiles
ggplot(cpue, aes(density)) +
geom_histogram() +
facet_wrap(~species2, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
cpue <- cpue %>%
group_by(species2) %>%
mutate(dens_quants = quantile(density, probs = 0.999)) %>%
filter(density < dens_quants) %>%
ungroup() %>%
dplyr::select(-dens_quants)
#> group_by: one grouping variable (species2)
#> mutate (grouped): new variable 'dens_quants' (double) with 3 unique values and 0% NA
#> filter (grouped): removed 30 rows (<1%), 27,867 rows remaining
#> ungroup: no grouping variables
# Make wide
wcpue <- cpue %>% pivot_wider(names_from = species2, values_from = density)
#> pivot_wider: reorganized (species2, density) into (flounder, large_cod, small_cod) [was 27867x3, now 9373x4]
head(wcpue)
#> # A tibble: 6 × 4
#> haul_id flounder large_cod small_cod
#> <chr> <dbl> <dbl> <dbl>
#> 1 1993:1:GFR:SOL:H20:21:1 9.41 77.9 0.628
#> 2 1993:1:GFR:SOL:H20:22:32 6.63 5.13 0
#> 3 1993:1:GFR:SOL:H20:23:31 0.953 2.61 0.00459
#> 4 1993:1:GFR:SOL:H20:24:30 1.89 9.71 0
#> 5 1993:1:GFR:SOL:H20:25:2 16.7 400. 4.78
#> 6 1993:1:GFR:SOL:H20:26:3 13.5 179. 2.95
# Left join in trawl information
nrow(wcpue)
#> [1] 9373
hauls <- cpue_full %>% distinct(haul_id, .keep_all = TRUE) %>% dplyr::select(-density)
#> distinct: removed 1,750,236 rows (99%), 9,373 rows remaining
nrow(hauls)
#> [1] 9373
cpue2 <- left_join(wcpue, hauls)
#> Joining, by = "haul_id"
#> left_join: added 20 columns (year, lat, lon, quarter, country, …)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 9,373
#> > =======
#> > rows total 9,373
colnames(cpue2)
#> [1] "haul_id" "flounder" "large_cod" "small_cod"
#> [5] "year" "lat" "lon" "quarter"
#> [9] "country" "month" "i_dx" "ices_rect"
#> [13] "sub_div" "length_cm" "id_haul_stomach" "species"
#> [17] "haul_id_size" "substrate" "depth" "temp"
#> [21] "oxy" "sal" "X" "Y"
cpue2 <- cpue2 %>% drop_na(oxy, temp, depth, sal, substrate, flounder, small_cod, large_cod)
#> drop_na: removed 405 rows (4%), 8,968 rows remaining
# Scale variables
cpue2 <- cpue2 %>%
mutate(depth_ct = depth - mean(depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(depth_sq)) / sd(depth_sq),
depth_sc = (depth - mean(depth)) / sd(depth),
temp_ct = temp - mean(temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(temp_sq)) / sd(temp_sq),
temp_sc = (temp - mean(temp)) / sd(temp),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
sal_sc = (sal - mean(sal)) / sd(sal),
flounder_sc = (flounder - mean(flounder) / sd(flounder)),
large_cod_sc = (large_cod - mean(large_cod) / sd(large_cod)),
small_cod_sc = (small_cod - mean(small_cod) / sd(small_cod)),
fyear = as.factor(year),
fsubstrate = as.factor(substrate))
#> mutate: new variable 'depth_ct' (double) with 149 unique values and 0% NA
#> new variable 'depth_sq' (double) with 149 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 149 unique values and 0% NA
#> new variable 'depth_sc' (double) with 149 unique values and 0% NA
#> new variable 'temp_ct' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sq' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 8,679 unique values and 0% NA
#> new variable 'temp_sc' (double) with 8,679 unique values and 0% NA
#> new variable 'oxy_sc' (double) with 8,678 unique values and 0% NA
#> new variable 'sal_sc' (double) with 8,690 unique values and 0% NA
#> new variable 'flounder_sc' (double) with 7,933 unique values and 0% NA
#> new variable 'large_cod_sc' (double) with 7,855 unique values and 0% NA
#> new variable 'small_cod_sc' (double) with 6,617 unique values and 0% NA
#> new variable 'fyear' (factor) with 28 unique values and 0% NA
#> new variable 'fsubstrate' (factor) with 5 unique values and 0% NA# colnames(cpue2)
#
# ggplot(cpue2, aes(flounder, large_cod)) +
# geom_point()
#
# ggplot(cpue2, aes(flounder, small_cod)) +
# geom_point()
#
# write.csv(cpue2, "output/mich_catch_by_length_q1_q4.csv")Read and scale the prediction grid
pred_grid_1 <- read_csv("data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- bind_rows(pred_grid_1, pred_grid_2)
# Scale with respect to data!
pred_grid <- pred_grid %>%
drop_na(substrate) %>%
mutate(X = X / 1000,
Y = Y / 1000,
depth_ct = depth - mean(cpue2$depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
temp_ct = temp - mean(cpue2$temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
oxy_sc = (oxy - mean(cpue2$oxy)) / sd(cpue2$oxy),
sal_sc = (sal - mean(cpue2$sal)) / sd(cpue2$sal),
fyear = as.factor(year),
fsubstrate = as.factor(substrate))
#> drop_na: removed 280 rows (<1%), 592,760 rows remaining
#> mutate: changed 592,760 values (100%) of 'X' (0 new NA)
#> changed 592,760 values (100%) of 'Y' (0 new NA)
#> new variable 'depth_ct' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sq' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 7,528 unique values and 0% NA
#> new variable 'depth_sc' (double) with 7,528 unique values and 0% NA
#> new variable 'temp_ct' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sq' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'temp_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'oxy_sc' (double) with 592,704 unique values and 0% NA
#> new variable 'sal_sc' (double) with 592,760 unique values and 0% NA
#> new variable 'fyear' (factor) with 28 unique values and 0% NA
#> new variable 'fsubstrate' (factor) with 5 unique values and 0% NA
# Split by quarter
pred_grid_q1 <- pred_grid %>% filter(quarter == 1)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
pred_grid_q4 <- pred_grid %>% filter(quarter == 4)
#> filter: removed 296,380 rows (50%), 296,380 rows remaining# Load models
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/spatial_overlap_cache/html")Plot data in space
# hist(cpue2$depth)
# hist(log(cpue2$depth))
#
# cpue_long <- cpue2 %>%
# pivot_longer(c(flounder, small_cod, large_cod)) %>%
# rename(species2 = name, density = value) %>%
# pivot_longer(c(depth, temp, oxy, sal)) %>%
# rename(env_var = name, env_var_value = value)
#
# ggplot(cpue_long, aes(density)) +
# geom_histogram() +
# facet_wrap(~species2, scales = "free", ncol = 1)
#
# ggplot(cpue_long, aes(env_var_value, density)) +
# geom_point(alpha = 0.3) +
# stat_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) +
# stat_smooth(method = "lm", formula = y ~ x + I(x^2), color = "red", se = FALSE) +
# facet_wrap(env_var~species2, scales = "free", ncol = 3) +
# coord_cartesian(expand = 0)
#
# cpue2 <- cpue2 %>% mutate(fle_presence = ifelse(flounder == 0, "N", "Y"),
# l_cod_presence = ifelse(large_cod == 0, "N", "Y"),
# s_cod_presence = ifelse(small_cod == 0, "N", "Y"))
#
# # Biomass density
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = flounder)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = large_cod)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = small_cod)) +
# scale_color_viridis_c(trans = "sqrt") +
# facet_wrap(~year)
#
# # Presence / absence
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = fle_presence)) +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = l_cod_presence)) +
# facet_wrap(~year)
#
# plot_map_fc +
# geom_point(data = cpue2, aes(X*1000, Y*1000, color = s_cod_presence)) +
# facet_wrap(~year)Create meshes
# Both Q's
spde <- make_mesh(cpue2,
xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans",
seed = 42)
# Q1
spde_q1 <- make_mesh(filter(cpue2, quarter == 1),
xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans",
seed = 42)
#> filter: removed 3,638 rows (41%), 5,330 rows remaining
# Q4
spde_q4 <- make_mesh(filter(cpue2, quarter == 4),
xy_cols = c("X", "Y"),
n_knots = 150,
type = "kmeans",
seed = 42)
#> filter: removed 5,330 rows (59%), 3,638 rows remainingFit models
q1
# Small cod model q1 spatial
mcod_s_q1_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining
sanity(mcod_s_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock -0.54972097 1.08136018 -2.6691480 1.5697060
#> 2 fsubstratehard clay 0.17092271 0.67999722 -1.1618474 1.5036928
#> 3 fsubstratehard-bottom complex 0.03007446 0.70107772 -1.3440126 1.4041615
#> 4 fsubstratemud -0.06888319 0.68011532 -1.4018847 1.2641183
#> 5 fsubstratesand -0.18812057 0.68131988 -1.5234830 1.1472418
#> 6 fyear1994 -0.27083562 0.65699125 -1.5585148 1.0168436
#> 7 fyear1995 0.32913371 0.64556718 -0.9361547 1.5944221
#> 8 fyear1996 0.20377555 0.65459642 -1.0792099 1.4867610
#> 9 fyear1997 -0.78333814 0.65633408 -2.0697293 0.5030530
#> 10 fyear1998 -0.15049380 0.64312981 -1.4110051 1.1100175
#> 11 fyear1999 -0.34980464 0.63759208 -1.5994622 0.8998529
#> 12 fyear2000 0.68322725 0.66498092 -0.6201114 1.9865659
#> 13 fyear2001 0.53330300 0.63057208 -0.7025956 1.7692016
#> 14 fyear2002 1.87881965 0.63487193 0.6344935 3.1231458
#> 15 fyear2003 0.25683392 0.65287546 -1.0227785 1.5364463
#> 16 fyear2004 -0.32113630 0.62569530 -1.5474766 0.9052040
#> 17 fyear2005 1.94428972 0.60379069 0.7608817 3.1276977
#> 18 fyear2006 0.17607793 0.61892330 -1.0369895 1.3891453
#> 19 fyear2007 0.88916597 0.60511521 -0.2968381 2.0751700
#> 20 fyear2008 1.22420903 0.60332687 0.0417101 2.4067080
#> 21 fyear2009 0.74141683 0.60380274 -0.4420148 1.9248485
#> 22 fyear2010 1.42101550 0.61003687 0.2253652 2.6166658
#> 23 fyear2011 0.23229999 0.62086565 -0.9845743 1.4491743
#> 24 fyear2012 0.33796827 0.61645874 -0.8702687 1.5462052
#> 25 fyear2013 1.49256114 0.60443447 0.3078914 2.6772309
#> 26 fyear2014 1.00617456 0.61271831 -0.1947313 2.2070804
#> 27 fyear2015 -0.24874693 0.61446510 -1.4530764 0.9555825
#> 28 fyear2016 -0.15828486 0.61471050 -1.3630953 1.0465256
#> 29 fyear2017 -0.30106379 0.61950192 -1.5152652 0.9131377
#> 30 fyear2018 0.64103439 0.61436656 -0.5631019 1.8451707
#> 31 fyear2019 -0.15021962 0.66175918 -1.4472438 1.1468045
#> 32 fyear2020 -1.20443537 0.66185914 -2.5016554 0.0927847
#> 33 depth_sc 0.97037576 0.08787148 0.7981508 1.1426007
#> 34 depth_sq_sc -1.28394502 0.07550417 -1.4319305 -1.1359596
#> 35 temp_sc 1.01359298 0.13782015 0.7434705 1.2837155
#> 36 temp_sq_sc -0.45935493 0.07285986 -0.6021576 -0.3165522
#> 37 oxy_sc 0.58282738 0.08955571 0.4073014 0.7583533
#> 38 sal_sc -0.08035481 0.10958924 -0.2951458 0.1344362# Large cod model q1 spatial
mcod_l_q1_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining
sanity(mcod_l_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 3.33190554 0.76118613 1.84000815 4.82380294
#> 2 fsubstratehard clay 3.96973799 0.47458439 3.03956967 4.89990631
#> 3 fsubstratehard-bottom complex 3.72985871 0.48918284 2.77107796 4.68863947
#> 4 fsubstratemud 3.96832494 0.47482006 3.03769473 4.89895516
#> 5 fsubstratesand 3.98132186 0.47566185 3.04904177 4.91360196
#> 6 fyear1994 0.40602594 0.52096627 -0.61504919 1.42710107
#> 7 fyear1995 0.12762982 0.51960316 -0.89077366 1.14603331
#> 8 fyear1996 0.97962265 0.52204960 -0.04357576 2.00282106
#> 9 fyear1997 -0.56368603 0.52591441 -1.59445932 0.46708727
#> 10 fyear1998 -0.64602778 0.52346347 -1.67199734 0.37994177
#> 11 fyear1999 -0.80355238 0.51875960 -1.82030251 0.21319775
#> 12 fyear2000 -0.46274022 0.54747744 -1.53577629 0.61029585
#> 13 fyear2001 0.04162165 0.51879413 -0.97519617 1.05843947
#> 14 fyear2002 0.28583472 0.52835100 -0.74971420 1.32138365
#> 15 fyear2003 0.26971748 0.53419605 -0.77728754 1.31672249
#> 16 fyear2004 -1.53828619 0.51955950 -2.55660410 -0.51996827
#> 17 fyear2005 -0.12409882 0.50031173 -1.10469180 0.85649416
#> 18 fyear2006 0.35259729 0.50247871 -0.63224289 1.33743746
#> 19 fyear2007 -0.28815905 0.50049926 -1.26911957 0.69280147
#> 20 fyear2008 0.11722341 0.49765679 -0.85816597 1.09261280
#> 21 fyear2009 0.38764681 0.49614749 -0.58478440 1.36007803
#> 22 fyear2010 0.34810389 0.50198280 -0.63576431 1.33197210
#> 23 fyear2011 0.09516022 0.50543903 -0.89548208 1.08580252
#> 24 fyear2012 -0.45146797 0.50800862 -1.44714657 0.54421063
#> 25 fyear2013 -0.13720362 0.49890510 -1.11503965 0.84063241
#> 26 fyear2014 -0.23631283 0.50478973 -1.22568251 0.75305686
#> 27 fyear2015 -0.16122928 0.50261088 -1.14632850 0.82386995
#> 28 fyear2016 -0.22682356 0.50268367 -1.21206545 0.75841833
#> 29 fyear2017 -0.54635722 0.50294671 -1.53211465 0.43940021
#> 30 fyear2018 -0.38825368 0.50554395 -1.37910160 0.60259425
#> 31 fyear2019 -1.02858635 0.54573504 -2.09820738 0.04103467
#> 32 fyear2020 -1.30537254 0.53538228 -2.35470252 -0.25604256
#> 33 depth_sc 0.75065191 0.06419549 0.62483106 0.87647276
#> 34 depth_sq_sc -0.39245921 0.03479079 -0.46064789 -0.32427052
#> 35 temp_sc 0.70706082 0.09830480 0.51438696 0.89973468
#> 36 temp_sq_sc -0.14525828 0.05124781 -0.24570215 -0.04481441
#> 37 oxy_sc 0.29743177 0.06620641 0.16766958 0.42719395
#> 38 sal_sc -0.17440207 0.07928522 -0.32979825 -0.01900590# Flounder model q1 spatial
mfle_q1_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 1),
mesh = spde_q1,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 3,638 rows (41%), 5,330 rows remaining
sanity(mfle_q1_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q1_space, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 4.09133618 0.52049350 3.071187672
#> 2 fsubstratehard clay 3.70825073 0.43066415 2.864164504
#> 3 fsubstratehard-bottom complex 3.60123785 0.43951943 2.739795601
#> 4 fsubstratemud 3.78796875 0.43054351 2.944118973
#> 5 fsubstratesand 3.68742673 0.43151917 2.841664697
#> 6 fyear1994 -0.10648718 0.41370545 -0.917334955
#> 7 fyear1995 0.03001146 0.41141847 -0.776353922
#> 8 fyear1996 0.07201784 0.41606651 -0.743457541
#> 9 fyear1997 -0.63738865 0.41839563 -1.457429017
#> 10 fyear1998 -0.74615179 0.41391762 -1.557415421
#> 11 fyear1999 -0.43988232 0.40810908 -1.239761423
#> 12 fyear2000 -0.03650738 0.43029747 -0.879874913
#> 13 fyear2001 0.07987207 0.41206573 -0.727761925
#> 14 fyear2002 0.37427374 0.41743409 -0.443882046
#> 15 fyear2003 0.28430880 0.41847879 -0.535894551
#> 16 fyear2004 -0.95578850 0.40545811 -1.750471793
#> 17 fyear2005 0.05309166 0.39438861 -0.719895802
#> 18 fyear2006 0.60071762 0.39551691 -0.174481285
#> 19 fyear2007 0.34827399 0.39233210 -0.420682810
#> 20 fyear2008 0.27884936 0.39317426 -0.491758030
#> 21 fyear2009 0.14276616 0.39269484 -0.626901588
#> 22 fyear2010 0.66332645 0.39474908 -0.110367540
#> 23 fyear2011 0.21630668 0.39770089 -0.563172740
#> 24 fyear2012 0.10076256 0.39578708 -0.674965867
#> 25 fyear2013 0.41927358 0.39221728 -0.349458156
#> 26 fyear2014 0.19635692 0.39694484 -0.581640683
#> 27 fyear2015 -0.27416323 0.39674972 -1.051778392
#> 28 fyear2016 -0.12358694 0.39610329 -0.899935112
#> 29 fyear2017 -0.08079854 0.39553248 -0.856027961
#> 30 fyear2018 0.26841628 0.39925996 -0.514118857
#> 31 fyear2019 -0.14305015 0.42921210 -0.984290404
#> 32 fyear2020 -0.29538565 0.41996741 -1.118506654
#> 33 depth_sc 0.09669746 0.05062436 -0.002524471
#> 34 depth_sq_sc -0.09632404 0.02637195 -0.148012125
#> 35 temp_sc 0.81118129 0.08580203 0.643012399
#> 36 temp_sq_sc -0.02569753 0.04448537 -0.112887244
#> 37 oxy_sc 0.51694645 0.05389546 0.411313283
#> 38 sal_sc 0.25730289 0.07541040 0.109501221
#> conf.high
#> 1 5.11148469
#> 2 4.55233697
#> 3 4.46268010
#> 4 4.63181852
#> 5 4.53318877
#> 6 0.70436060
#> 7 0.83637685
#> 8 0.88749322
#> 9 0.18265172
#> 10 0.06511185
#> 11 0.35999679
#> 12 0.80686016
#> 13 0.88750606
#> 14 1.19242952
#> 15 1.10451215
#> 16 -0.16110521
#> 17 0.82607913
#> 18 1.37591652
#> 19 1.11723078
#> 20 1.04945675
#> 21 0.91243391
#> 22 1.43702043
#> 23 0.99578610
#> 24 0.87649099
#> 25 1.18800532
#> 26 0.97435451
#> 27 0.50345194
#> 28 0.65276124
#> 29 0.69443088
#> 30 1.05095142
#> 31 0.69819011
#> 32 0.52773536
#> 33 0.19591939
#> 34 -0.04463596
#> 35 0.97935018
#> 36 0.06149219
#> 37 0.62257961
#> 38 0.40510457q4
# Small cod model q4 spatial
mcod_s_q4_space <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining
sanity(mcod_s_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 0.32249199 0.73538955 -1.1188450 1.76382903
#> 2 fsubstratehard clay -0.62378219 0.55544819 -1.7124406 0.46487626
#> 3 fsubstratehard-bottom complex -0.96634332 0.57126621 -2.0860045 0.15331788
#> 4 fsubstratemud -0.89407109 0.55349970 -1.9789106 0.19076839
#> 5 fsubstratesand -0.67373306 0.55549607 -1.7624853 0.41501922
#> 6 fyear1994 0.06757994 0.68460289 -1.2742171 1.40937695
#> 7 fyear1995 1.05387153 0.65331032 -0.2265932 2.33433623
#> 8 fyear1996 -0.74286053 0.69564320 -2.1062961 0.62057508
#> 9 fyear1997 0.54977163 0.63147742 -0.6879014 1.78744463
#> 10 fyear1998 0.31989404 0.65061515 -0.9552882 1.59507631
#> 11 fyear1999 0.51812992 0.63184376 -0.7202611 1.75652093
#> 12 fyear2000 1.07924843 0.60720418 -0.1108499 2.26934676
#> 13 fyear2001 2.04473353 0.59447261 0.8795886 3.20987843
#> 14 fyear2002 1.44874905 0.59584509 0.2809141 2.61658397
#> 15 fyear2003 -0.30863035 0.61053166 -1.5052504 0.88798971
#> 16 fyear2004 2.37016410 0.58870662 1.2163203 3.52400788
#> 17 fyear2005 1.47212495 0.57301842 0.3490295 2.59522040
#> 18 fyear2006 1.83767110 0.56994858 0.7205924 2.95474978
#> 19 fyear2007 2.00583510 0.56804833 0.8924808 3.11918938
#> 20 fyear2008 1.54410451 0.56548803 0.4357683 2.65244068
#> 21 fyear2009 2.26931005 0.56514325 1.1616496 3.37697046
#> 22 fyear2010 1.57705454 0.56759490 0.4645890 2.68952010
#> 23 fyear2011 0.75334444 0.57103716 -0.3658678 1.87255671
#> 24 fyear2012 1.31095911 0.57701625 0.1800280 2.44189018
#> 25 fyear2013 1.86685633 0.57015276 0.7493775 2.98433522
#> 26 fyear2014 0.73565198 0.57489336 -0.3911183 1.86242226
#> 27 fyear2015 0.33706211 0.57855532 -0.7968855 1.47100970
#> 28 fyear2016 -0.01423134 0.57677617 -1.1446919 1.11622918
#> 29 fyear2017 1.41721173 0.56912509 0.3017471 2.53267640
#> 30 fyear2018 0.37104157 0.57473009 -0.7554087 1.49749184
#> 31 fyear2019 -0.15607354 0.63784335 -1.4062235 1.09407646
#> 32 fyear2020 0.12153068 0.60740905 -1.0689692 1.31203054
#> 33 depth_sc 0.28276560 0.09106935 0.1042730 0.46125823
#> 34 depth_sq_sc -1.36146413 0.08572625 -1.5294845 -1.19344378
#> 35 temp_sc 0.82467933 0.13390142 0.5622374 1.08712129
#> 36 temp_sq_sc -0.20336686 0.07336668 -0.3471629 -0.05957081
#> 37 oxy_sc 0.99602212 0.10423533 0.7917246 1.20031962
#> 38 sal_sc 0.46967220 0.11623198 0.2418617 0.69748268# Large cod model q4 spatial
mcod_l_q4_space <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining
sanity(mcod_l_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 4.461261842 0.67727269 3.13383176
#> 2 fsubstratehard clay 4.085728133 0.50076507 3.10424663
#> 3 fsubstratehard-bottom complex 3.766673734 0.51197346 2.76322418
#> 4 fsubstratemud 3.913855555 0.49969560 2.93447018
#> 5 fsubstratesand 3.887435320 0.50143270 2.90464528
#> 6 fyear1994 0.419584087 0.56203685 -0.68198790
#> 7 fyear1995 1.069600130 0.55806448 -0.02418614
#> 8 fyear1996 0.257708655 0.56677971 -0.85315916
#> 9 fyear1997 -0.370530757 0.54713854 -1.44290259
#> 10 fyear1998 -0.204480195 0.55538387 -1.29301257
#> 11 fyear1999 -0.631959902 0.55349467 -1.71678952
#> 12 fyear2000 0.073268003 0.53256216 -0.97053465
#> 13 fyear2001 0.160909371 0.52872110 -0.87536495
#> 14 fyear2002 0.776504998 0.52166877 -0.24594700
#> 15 fyear2003 -1.483437852 0.53590092 -2.53378435
#> 16 fyear2004 0.096387087 0.52463275 -0.93187420
#> 17 fyear2005 0.600156755 0.50644354 -0.39245434
#> 18 fyear2006 0.332251050 0.50594943 -0.65939161
#> 19 fyear2007 0.512317962 0.50465526 -0.47678817
#> 20 fyear2008 0.655240592 0.50146897 -0.32762053
#> 21 fyear2009 0.513073198 0.50309145 -0.47296793
#> 22 fyear2010 0.676087205 0.50348831 -0.31073176
#> 23 fyear2011 -0.007275424 0.50555201 -0.99813915
#> 24 fyear2012 -0.396814061 0.51432782 -1.40487806
#> 25 fyear2013 0.002090184 0.50732834 -0.99225509
#> 26 fyear2014 -0.202627548 0.50870179 -1.19966474
#> 27 fyear2015 -0.087591060 0.50823553 -1.08371439
#> 28 fyear2016 -0.818491405 0.51022386 -1.81851179
#> 29 fyear2017 -0.424939151 0.50767462 -1.41996313
#> 30 fyear2018 -1.153193343 0.51144123 -2.15559973
#> 31 fyear2019 -1.254966349 0.55782243 -2.34827821
#> 32 fyear2020 -1.536938007 0.54214527 -2.59952322
#> 33 depth_sc 0.408087994 0.07078593 0.26935012
#> 34 depth_sq_sc -0.542902793 0.05660972 -0.65385581
#> 35 temp_sc 0.239683585 0.10150667 0.04073417
#> 36 temp_sq_sc -0.063530409 0.05595097 -0.17319229
#> 37 oxy_sc 0.722704596 0.07857480 0.56870083
#> 38 sal_sc 0.288968482 0.09016190 0.11225441
#> conf.high
#> 1 5.78869192
#> 2 5.06720964
#> 3 4.77012329
#> 4 4.89324093
#> 5 4.87022536
#> 6 1.52115608
#> 7 2.16338640
#> 8 1.36857647
#> 9 0.70184108
#> 10 0.88405218
#> 11 0.45286972
#> 12 1.11707066
#> 13 1.19718369
#> 14 1.79895699
#> 15 -0.43309135
#> 16 1.12464838
#> 17 1.59276785
#> 18 1.32389371
#> 19 1.50142409
#> 20 1.63810171
#> 21 1.49911433
#> 22 1.66290617
#> 23 0.98358831
#> 24 0.61124994
#> 25 0.99643546
#> 26 0.79440965
#> 27 0.90853227
#> 28 0.18152898
#> 29 0.57008483
#> 30 -0.15078695
#> 31 -0.16165448
#> 32 -0.47435279
#> 33 0.54682587
#> 34 -0.43194978
#> 35 0.43863300
#> 36 0.04613147
#> 37 0.87670836
#> 38 0.46568256# Flounder model q4 spatial
mfle_q4_space <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = filter(cpue2, quarter == 4),
mesh = spde_q4,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> filter: removed 5,330 rows (59%), 3,638 rows remaining
sanity(mfle_q4_space)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle_q4_space, conf.int = TRUE)
#> term estimate std.error conf.low conf.high
#> 1 fsubstratebedrock 2.71204220 0.60223670 1.53167996 3.89240444
#> 2 fsubstratehard clay 2.09065120 0.56330160 0.98660035 3.19470205
#> 3 fsubstratehard-bottom complex 1.98083675 0.56839418 0.86680462 3.09486888
#> 4 fsubstratemud 2.01546063 0.56224288 0.91348483 3.11743643
#> 5 fsubstratesand 2.10477305 0.56324257 1.00083789 3.20870820
#> 6 fyear1994 -0.42238110 0.59608389 -1.59068405 0.74592184
#> 7 fyear1995 0.82691107 0.58064409 -0.31113043 1.96495257
#> 8 fyear1996 0.34292251 0.59481583 -0.82289509 1.50874011
#> 9 fyear1997 0.24815800 0.56532853 -0.85986556 1.35618156
#> 10 fyear1998 -0.26640178 0.57547313 -1.39430839 0.86150484
#> 11 fyear1999 -0.96023820 0.57740366 -2.09192858 0.17145219
#> 12 fyear2000 0.06955067 0.55504887 -1.01832512 1.15742646
#> 13 fyear2001 0.14979012 0.55147291 -0.93107693 1.23065717
#> 14 fyear2002 0.87163543 0.54406255 -0.19470757 1.93797844
#> 15 fyear2003 -0.99038344 0.54565467 -2.05984694 0.07908007
#> 16 fyear2004 0.51169827 0.54311299 -0.55278363 1.57618016
#> 17 fyear2005 0.52266141 0.52661866 -0.50949221 1.55481502
#> 18 fyear2006 0.58363248 0.52306129 -0.44154881 1.60881377
#> 19 fyear2007 0.68479501 0.52247358 -0.33923438 1.70882440
#> 20 fyear2008 0.88301325 0.51990783 -0.13598738 1.90201388
#> 21 fyear2009 0.90246533 0.51882379 -0.11441062 1.91934127
#> 22 fyear2010 1.21271007 0.51853304 0.19640400 2.22901615
#> 23 fyear2011 0.87184325 0.51905899 -0.14549368 1.88918019
#> 24 fyear2012 0.48386735 0.52304030 -0.54127280 1.50900751
#> 25 fyear2013 0.71014527 0.52389239 -0.31666494 1.73695548
#> 26 fyear2014 0.42636457 0.52325026 -0.59918710 1.45191624
#> 27 fyear2015 0.49211851 0.52179361 -0.53057818 1.51481520
#> 28 fyear2016 0.36022963 0.52073115 -0.66038467 1.38084393
#> 29 fyear2017 0.56235857 0.52240042 -0.46152744 1.58624458
#> 30 fyear2018 0.37786922 0.52380208 -0.64876400 1.40450244
#> 31 fyear2019 0.58021712 0.56227654 -0.52182464 1.68225888
#> 32 fyear2020 0.03493564 0.54492131 -1.03309051 1.10296178
#> 33 depth_sc -0.06312977 0.06548158 -0.19147132 0.06521178
#> 34 depth_sq_sc -0.83253692 0.06332174 -0.95664524 -0.70842859
#> 35 temp_sc 0.57487192 0.10083785 0.37723337 0.77251048
#> 36 temp_sq_sc -0.03971065 0.05665936 -0.15076096 0.07133966
#> 37 oxy_sc 1.20536646 0.08323146 1.04223581 1.36849712
#> 38 sal_sc 0.16223809 0.09127951 -0.01666646 0.34114264Models with the other species as coeffcients
We can’t put them in the main model because we want to use that for spatial overlap. There isn’t as much need to fit the models below by quarter (though we could), because while the distribution changes, we don’t know a priori that the effect should.
# Small cod model
mcod_s <- sdmTMB(small_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc + flounder_sc,
data = cpue2,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation
sanity(mcod_s)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_s, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 0.628988603 0.6284130899 -0.60267842
#> 2 fsubstratehard clay 0.186705195 0.4303576376 -0.65678028
#> 3 fsubstratehard-bottom complex 0.139329862 0.4463497732 -0.73549962
#> 4 fsubstratemud 0.021906096 0.4298307553 -0.82054670
#> 5 fsubstratesand 0.045164068 0.4307133542 -0.79901859
#> 6 fyear1994 -0.647472621 0.4659866755 -1.56078972
#> 7 fyear1995 0.135138774 0.4538600093 -0.75441050
#> 8 fyear1996 -0.585039168 0.4663404873 -1.49904973
#> 9 fyear1997 -0.781861845 0.4576482420 -1.67883592
#> 10 fyear1998 -0.456436267 0.4538632659 -1.34599192
#> 11 fyear1999 -0.475663483 0.4499229927 -1.35749634
#> 12 fyear2000 0.363354134 0.4566324637 -0.53162905
#> 13 fyear2001 0.581528036 0.4369675478 -0.27491262
#> 14 fyear2002 1.002107172 0.4341469713 0.15119474
#> 15 fyear2003 -0.503903733 0.4517292266 -1.38927675
#> 16 fyear2004 0.432451462 0.4312379975 -0.41275948
#> 17 fyear2005 1.149304640 0.4156855586 0.33457592
#> 18 fyear2006 0.496325702 0.4223847818 -0.33153326
#> 19 fyear2007 0.759332204 0.4169238409 -0.05782351
#> 20 fyear2008 0.771594631 0.4156141740 -0.04299418
#> 21 fyear2009 0.990532701 0.4140889987 0.17893318
#> 22 fyear2010 0.831369822 0.4184572443 0.01120869
#> 23 fyear2011 -0.040597968 0.4221376512 -0.86797256
#> 24 fyear2012 0.251791028 0.4241668674 -0.57956076
#> 25 fyear2013 1.225870027 0.4168090620 0.40893928
#> 26 fyear2014 0.546055609 0.4232892680 -0.28357611
#> 27 fyear2015 -0.426740925 0.4262169679 -1.26211083
#> 28 fyear2016 -0.370330338 0.4232011663 -1.19978938
#> 29 fyear2017 0.309418397 0.4208926310 -0.51551600
#> 30 fyear2018 0.135837741 0.4199787271 -0.68730544
#> 31 fyear2019 -0.404787529 0.4623097601 -1.31089801
#> 32 fyear2020 -1.015116571 0.4529117231 -1.90280724
#> 33 depth_sc 0.563637368 0.0671644517 0.43199746
#> 34 depth_sq_sc -1.086366763 0.0591121348 -1.20222442
#> 35 temp_sc 0.620151439 0.0304481960 0.56047407
#> 36 temp_sq_sc -0.380384042 0.0230597138 -0.42558025
#> 37 oxy_sc 0.506692384 0.0572793122 0.39442700
#> 38 sal_sc 0.070215044 0.0591226458 -0.04566321
#> 39 flounder_sc 0.002697162 0.0001541721 0.00239499
#> conf.high
#> 1 1.860655627
#> 2 1.030190665
#> 3 1.014159342
#> 4 0.864358896
#> 5 0.889346730
#> 6 0.265844480
#> 7 1.024688046
#> 8 0.328971392
#> 9 0.115112227
#> 10 0.433119388
#> 11 0.406169378
#> 12 1.258337317
#> 13 1.437968692
#> 14 1.853019599
#> 15 0.381469281
#> 16 1.277662405
#> 17 1.964033363
#> 18 1.324184662
#> 19 1.576487917
#> 20 1.586183444
#> 21 1.802132225
#> 22 1.651530950
#> 23 0.786776625
#> 24 1.083142811
#> 25 2.042800777
#> 26 1.375687329
#> 27 0.408628982
#> 28 0.459128706
#> 29 1.134352795
#> 30 0.958980920
#> 31 0.501322951
#> 32 -0.127425906
#> 33 0.695277274
#> 34 -0.970509107
#> 35 0.679828807
#> 36 -0.335187834
#> 37 0.618957773
#> 38 0.186093301
#> 39 0.002999333# Large cod model q4 spatial
mcod_l <- sdmTMB(large_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc + flounder_sc,
data = cpue2,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation
sanity(mcod_l)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mcod_l, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 3.42710819 0.5776955169 2.294845779
#> 2 fsubstratehard clay 3.65473107 0.4329649614 2.806135343
#> 3 fsubstratehard-bottom complex 3.53484686 0.4407693938 2.670954722
#> 4 fsubstratemud 3.63725350 0.4326850367 2.789206409
#> 5 fsubstratesand 3.61916586 0.4331881805 2.770132626
#> 6 fyear1994 0.19742496 0.4461982359 -0.677107517
#> 7 fyear1995 0.17722202 0.4452049165 -0.695363579
#> 8 fyear1996 0.59342116 0.4467822919 -0.282256041
#> 9 fyear1997 -0.62579353 0.4465154309 -1.500947688
#> 10 fyear1998 -0.58927473 0.4457547961 -1.462938074
#> 11 fyear1999 -0.78880418 0.4439619893 -1.658953688
#> 12 fyear2000 -0.16735879 0.4541355806 -1.057448169
#> 13 fyear2001 -0.10073693 0.4411441110 -0.965363495
#> 14 fyear2002 0.29062836 0.4386040328 -0.569019749
#> 15 fyear2003 -0.59834008 0.4507702144 -1.481833463
#> 16 fyear2004 -0.68491754 0.4355127979 -1.538506941
#> 17 fyear2005 0.16123839 0.4238126275 -0.669419100
#> 18 fyear2006 0.30196408 0.4258227266 -0.532633127
#> 19 fyear2007 -0.05063772 0.4254537613 -0.884511765
#> 20 fyear2008 0.23752649 0.4236463810 -0.592805160
#> 21 fyear2009 0.47468725 0.4219802066 -0.352378757
#> 22 fyear2010 0.47630118 0.4249369173 -0.356559877
#> 23 fyear2011 -0.10742847 0.4265686555 -0.943487677
#> 24 fyear2012 -0.45734127 0.4305572091 -1.301217898
#> 25 fyear2013 -0.06534330 0.4252864066 -0.898889341
#> 26 fyear2014 -0.10135916 0.4285620124 -0.941325273
#> 27 fyear2015 -0.18794203 0.4278915386 -1.026594032
#> 28 fyear2016 -0.55035878 0.4276565701 -1.388550257
#> 29 fyear2017 -0.52852784 0.4268257001 -1.365090842
#> 30 fyear2018 -0.72796029 0.4276077623 -1.566056107
#> 31 fyear2019 -1.15141866 0.4604934626 -2.053969258
#> 32 fyear2020 -1.47517542 0.4485627645 -2.354342285
#> 33 depth_sc 0.48942278 0.0492980160 0.392800447
#> 34 depth_sq_sc -0.31765342 0.0284809234 -0.373475003
#> 35 temp_sc 0.32629619 0.0213972972 0.284358257
#> 36 temp_sq_sc -0.25106483 0.0164175004 -0.283242543
#> 37 oxy_sc 0.23252038 0.0414262955 0.151326338
#> 38 sal_sc -0.14439376 0.0426861689 -0.228057118
#> 39 flounder_sc 0.00223198 0.0001122465 0.002011981
#> conf.high
#> 1 4.559370594
#> 2 4.503326805
#> 3 4.398738997
#> 4 4.485300587
#> 5 4.468199091
#> 6 1.071957428
#> 7 1.049807625
#> 8 1.469098361
#> 9 0.249360638
#> 10 0.284388618
#> 11 0.081345331
#> 12 0.722730595
#> 13 0.763889644
#> 14 1.150276466
#> 15 0.285153308
#> 16 0.168671857
#> 17 0.991895872
#> 18 1.136561289
#> 19 0.783236333
#> 20 1.067858138
#> 21 1.301753257
#> 22 1.309162231
#> 23 0.728630727
#> 24 0.386535348
#> 25 0.768202739
#> 26 0.738606946
#> 27 0.650709978
#> 28 0.287832693
#> 29 0.308035157
#> 30 0.110135520
#> 31 -0.248868055
#> 32 -0.596008558
#> 33 0.586045118
#> 34 -0.261831835
#> 35 0.368234121
#> 36 -0.218887124
#> 37 0.313714432
#> 38 -0.060730410
#> 39 0.002451979# Flounder model q4 spatial
mfle <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc + large_cod_sc + small_cod_sc,
data = cpue2,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "IID",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation
#> Warning in stats::nlminb(start = tmb_obj$par, objective = tmb_obj$fn, gradient =
#> tmb_obj$gr, : NA/NaN function evaluation
sanity(mfle)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfle, conf.int = TRUE)
#> term estimate std.error conf.low
#> 1 fsubstratebedrock 3.7480605290 4.809351e-01 2.8054450987
#> 2 fsubstratehard clay 3.3821101833 4.265386e-01 2.5461098645
#> 3 fsubstratehard-bottom complex 3.4792137737 4.328010e-01 2.6309394082
#> 4 fsubstratemud 3.5188278644 4.260961e-01 2.6836949097
#> 5 fsubstratesand 3.4790847848 4.272366e-01 2.6417164564
#> 6 fyear1994 -0.1746850881 3.764376e-01 -0.9124892955
#> 7 fyear1995 0.0322147643 3.746166e-01 -0.7020202906
#> 8 fyear1996 0.0789515893 3.779533e-01 -0.6618232722
#> 9 fyear1997 -0.5308645488 3.770456e-01 -1.2698602707
#> 10 fyear1998 -0.7059549558 3.762681e-01 -1.4434268281
#> 11 fyear1999 -0.5208658197 3.709189e-01 -1.2478534311
#> 12 fyear2000 -0.2732915952 3.851104e-01 -1.0280940453
#> 13 fyear2001 -0.0062696627 3.727832e-01 -0.7369112902
#> 14 fyear2002 0.2945432822 3.674595e-01 -0.4256641871
#> 15 fyear2003 -0.3504426411 3.734771e-01 -1.0824442460
#> 16 fyear2004 -0.3602718594 3.635209e-01 -1.0727596604
#> 17 fyear2005 -0.0703193795 3.544560e-01 -0.7650404426
#> 18 fyear2006 0.3092207848 3.551067e-01 -0.3867755284
#> 19 fyear2007 0.1500276245 3.534010e-01 -0.5426256724
#> 20 fyear2008 0.2362681297 3.535750e-01 -0.4567262252
#> 21 fyear2009 0.1724025532 3.531042e-01 -0.5196688931
#> 22 fyear2010 0.5688610670 3.539224e-01 -0.1248140049
#> 23 fyear2011 0.1819444743 3.544701e-01 -0.5128041571
#> 24 fyear2012 0.0013916538 3.562480e-01 -0.6968415743
#> 25 fyear2013 0.3524349089 3.537225e-01 -0.3408483580
#> 26 fyear2014 -0.1340147825 3.581498e-01 -0.8359754149
#> 27 fyear2015 -0.2295553924 3.567940e-01 -0.9288587597
#> 28 fyear2016 -0.2423025107 3.555494e-01 -0.9391664435
#> 29 fyear2017 -0.1301990053 3.548341e-01 -0.8256611383
#> 30 fyear2018 0.0039169536 3.565250e-01 -0.6948592487
#> 31 fyear2019 0.0590546671 3.854075e-01 -0.6963301850
#> 32 fyear2020 -0.4085526520 3.748436e-01 -1.1432325475
#> 33 depth_sc 0.1779467626 4.925036e-02 0.0814178360
#> 34 depth_sq_sc -0.0236437979 2.638278e-02 -0.0753531027
#> 35 temp_sc 0.4247440518 2.082157e-02 0.3839345247
#> 36 temp_sq_sc -0.3098803388 1.650299e-02 -0.3422256124
#> 37 oxy_sc 0.3911327055 4.277491e-02 0.3072954156
#> 38 sal_sc 0.2797308385 4.675097e-02 0.1881006181
#> 39 large_cod_sc 0.0005973076 4.037953e-05 0.0005181651
#> 40 small_cod_sc 0.0038476870 4.201599e-04 0.0030241887
#> conf.high
#> 1 4.690675959
#> 2 4.218110502
#> 3 4.327488139
#> 4 4.353960819
#> 5 4.316453113
#> 6 0.563119119
#> 7 0.766449819
#> 8 0.819726451
#> 9 0.208131173
#> 10 0.031516917
#> 11 0.206121792
#> 12 0.481510855
#> 13 0.724371965
#> 14 1.014750751
#> 15 0.381558964
#> 16 0.352215942
#> 17 0.624401684
#> 18 1.005217098
#> 19 0.842680922
#> 20 0.929262485
#> 21 0.864474000
#> 22 1.262536139
#> 23 0.876693106
#> 24 0.699624882
#> 25 1.045718176
#> 26 0.567945850
#> 27 0.469747975
#> 28 0.454561422
#> 29 0.565263128
#> 30 0.702693156
#> 31 0.814439519
#> 32 0.326127243
#> 33 0.274475689
#> 34 0.028065507
#> 35 0.465553579
#> 36 -0.277535065
#> 37 0.474969995
#> 38 0.371361059
#> 39 0.000676450
#> 40 0.004671185Plot conditional density effects
mfle_d <- visreg(mfle, xvar = "small_cod_sc", plot = FALSE)
mfle_d2 <- visreg(mfle, xvar = "large_cod_sc", plot = FALSE)
mscod_d <- visreg(mcod_s, xvar = "flounder_sc", plot = FALSE)
mlcod_d <- visreg(mcod_l, xvar = "flounder_sc", plot = FALSE)pal <- brewer.pal(n = 3, name = "Set1")[2]
p1 <- ggplot(mfle_d$fit, aes(x = small_cod_sc, y = visregFit)) +
geom_point(aes(y = visregRes), data = mfle_d$res, size = 1, alpha = 0.2, color = "grey30") +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
geom_line(color = pal, size = 1.2) +
labs(x = "Small cod density",
y = "Flounder visregFit")
p2 <- ggplot(mfle_d2$fit, aes(x = large_cod_sc, y = visregFit)) +
geom_point(aes(y = visregRes), data = mfle_d2$res, size = 1, alpha = 0.2, color = "grey30") +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
geom_line(color = pal, size = 1.2) +
labs(x = "Large cod density",
y = "Flounder visregFit")
p3 <- ggplot(mscod_d$fit, aes(x = flounder_sc, y = visregFit)) +
geom_point(aes(y = visregRes), data = mscod_d$res, size = 1, alpha = 0.2, color = "grey30") +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
geom_line(color = pal, size = 1.2) +
labs(x = "Flounder density",
y = "Small cod visregFit")
p4 <- ggplot(mlcod_d$fit, aes(x = flounder_sc, y = visregFit)) +
geom_point(aes(y = visregRes), data = mlcod_d$res, size = 1, alpha = 0.2, color = "grey30") +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2) +
geom_line(color = pal, size = 1.2) +
labs(x = "Flounder density",
y = "Large cod visregFit")
p1 + p2 + p3 + p4 & theme(aspect.ratio = 1)
ggsave("figures/density_coefs_smooth.pdf", width = 20, height = 20, units = "cm")Plot coefficients
# Q1
# small cod
small_cod_q1 <- tidy(mcod_s_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Small cod",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# large cod
large_cod_q1 <- tidy(mcod_l_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Large cod",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# flounder
flounder_q1 <- tidy(mfle_q1_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Flounder",
quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
# small cod
small_cod_q4 <- tidy(mcod_s_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Small cod",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# large cod
large_cod_q4 <- tidy(mcod_l_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Large cod",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# flounder
flounder_q4 <- tidy(mfle_q4_space, effects = "fixed", conf.int = TRUE) %>%
mutate(species = "Flounder",
quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
coefs <- bind_rows(small_cod_q1,
large_cod_q1,
flounder_q1,
small_cod_q4,
large_cod_q4,
flounder_q4) %>%
mutate(param_group = ifelse(term %in% c("fsubstratebedrock",
"fsubstratehard clay",
"fsubstratehard-bottom complex",
"fsubstratemud",
"fsubstratesand"),
"Substrate",
"Continious")) %>%
mutate(term = ifelse(term == "fsubstratebedrock", "Bedrock", term),
term = ifelse(term == "fsubstratehard clay", "Hard clay", term),
term = ifelse(term == "fsubstratehard-bottom complex", "Hard-bottom\ncomplex", term),
term = ifelse(term == "fsubstratemud", "Mud", term),
term = ifelse(term == "fsubstratesand", "Sand", term),
term = ifelse(term == "depth_sc", "Depth", term),
term = ifelse(term == "depth_sq_sc", "Depth^2", term),
term = ifelse(term == "temp_sc", "Temperature", term),
term = ifelse(term == "temp_sq_sc", "Temperature^2", term),
term = ifelse(term == "oxy_sc", "Oxygen", term),
term = ifelse(term == "sal_sc", "Salinity", term))
#> mutate: new variable 'param_group' (character) with 2 unique values and 0% NA
#> mutate: changed 66 values (29%) of 'term' (0 new NA)
coefs %>%
filter(!grepl('year', term)) %>%
ggplot(aes(term, estimate, color = species, shape = factor(quarter))) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 2) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
labs(shape = "Quarter") +
facet_wrap(~param_group, scales = "free") +
labs(x = "Predictor", y = "Standardized coefficient") +
theme_plot() +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
aspect.ratio = 1) +
NULL
#> filter: removed 162 rows (71%), 66 rows remaining
ggsave("figures/habitat_coefs.pdf", width = 20, height = 20, units = "cm")Predict on grid
mcod_s_q1_pred <- predict(mcod_s_q1_space, newdata = pred_grid_q1)
mcod_s_q4_pred <- predict(mcod_s_q4_space, newdata = pred_grid_q4)
mcod_l_q1_pred <- predict(mcod_l_q1_space, newdata = pred_grid_q1)
mcod_l_q4_pred <- predict(mcod_l_q4_space, newdata = pred_grid_q4)
mfle_q1_pred <- predict(mfle_q1_space, newdata = pred_grid_q1)
mfle_q4_pred <- predict(mfle_q4_space, newdata = pred_grid_q4)Plot predictions on map
p1 <- plot_map +
geom_raster(data = bind_rows(mcod_s_q1_pred, mcod_s_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Cod < 25 cm") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
p2 <- plot_map +
geom_raster(data = bind_rows(mcod_l_q1_pred, mcod_l_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Cod > 25 cm") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
p3 <- plot_map +
geom_raster(data = bind_rows(mfle_q1_pred, mfle_q4_pred) %>%
filter(year %in% c("1995", "2019")),
aes(X*1000, Y*1000, fill = exp(est))) +
facet_grid(quarter~year) +
ggtitle("Flounder") +
scale_fill_viridis(trans = "sqrt")
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
((p1 + p2) / (p3 + plot_spacer())) + plot_annotation(tag_levels = "A") &
theme(legend.text = element_text(angle = 90))
# Check flounder...
plot_map +
geom_raster(data = mfle_q1_pred,
aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt")
plot_map +
geom_raster(data = mfle_q4_pred,
aes(X*1000, Y*1000, fill = exp(est))) +
facet_wrap(~year) +
scale_fill_viridis(trans = "sqrt")
Calculate spatial overlap indices
all_pred_q1 <- mcod_s_q1_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA
all_pred_q1 <- all_pred_q1 %>%
mutate(large_cod = exp(mcod_l_q1_pred$est),
fle = exp(mfle_q1_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#> new variable 'fle' (double) with 296,380 unique values and 0% NA
all_pred_q4 <- mcod_s_q4_pred %>% mutate(small_cod = exp(est))
#> mutate: new variable 'small_cod' (double) with 296,380 unique values and 0% NA
all_pred_q4 <- all_pred_q4 %>%
mutate(large_cod = exp(mcod_l_q4_pred$est),
fle = exp(mfle_q4_pred$est))
#> mutate: new variable 'large_cod' (double) with 296,380 unique values and 0% NA
#> new variable 'fle' (double) with 296,380 unique values and 0% NA
all_pred <- bind_rows(all_pred_q1, all_pred_q4)
# Calculate biomass overlap
loc_collocspfn <- function(prey, pred) {
p_prey <- prey/sum(prey, na.rm = T)
p_pred <- pred/sum(pred, na.rm = T)
(p_prey*p_pred)/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}
loc_collocspfn_tot <- function(prey, pred) {
p_prey <- prey/sum(prey, na.rm = T)
p_pred <- pred/sum(pred, na.rm = T)
sum((p_prey*p_pred))/(sqrt(sum(p_prey^2, na.rm = T))* sqrt(sum(p_pred^2, na.rm = T)))
}
all_pred <- all_pred %>%
group_by(year, quarter) %>%
mutate(scod_fle_ovr = loc_collocspfn(pred = small_cod, prey = fle),
scod_fle_ovr_tot = loc_collocspfn_tot(pred = small_cod, prey = fle),
lcod_fle_ovr = loc_collocspfn(pred = large_cod, prey = fle),
lcod_fle_ovr_tot = loc_collocspfn_tot(pred = large_cod, prey = fle)) %>%
ungroup()
#> group_by: 2 grouping variables (year, quarter)
#> mutate (grouped): new variable 'scod_fle_ovr' (double) with 592,760 unique values and 0% NA
#> new variable 'scod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> new variable 'lcod_fle_ovr' (double) with 592,760 unique values and 0% NA
#> new variable 'lcod_fle_ovr_tot' (double) with 56 unique values and 0% NA
#> ungroup: no grouping variablesPlot
plot_map_fc +
geom_raster(data = all_pred %>% filter(quarter == 1),
aes(x = X, y = Y, fill = scod_fle_ovr)) +
scale_fill_viridis_c(trans = "sqrt", name = "small cod-flounder") +
facet_wrap(~ year, ncol = 5) +
theme(legend.text = element_text(angle = 90)) +
NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).
plot_map_fc +
geom_raster(data = all_pred %>% filter(quarter == 4),
aes(x = X, y = Y, fill = lcod_fle_ovr)) +
scale_fill_viridis_c(trans = "sqrt", name = "large cod-flounder") +
facet_wrap(~ year, ncol = 5) +
theme(legend.text = element_text(angle = 90)) +
NULL
#> filter: removed 296,380 rows (50%), 296,380 rows remaining
#> Warning: Removed 296380 rows containing missing values (geom_raster).
# In space for selected year
all_pred2 <- all_pred %>%
filter(year %in% c(1995, 2015)) %>%
dplyr::select(X, Y, scod_fle_ovr, lcod_fle_ovr, quarter, year) %>%
pivot_longer(c(scod_fle_ovr, lcod_fle_ovr)) %>%
rename(overlap = name) %>%
mutate(overlap = ifelse(overlap == "scod_fle_ovr", "Small cod : Flounder", overlap),
overlap = ifelse(overlap == "lcod_fle_ovr", "Large cod : Flounder", overlap))
#> filter: removed 550,420 rows (93%), 42,340 rows remaining
#> pivot_longer: reorganized (scod_fle_ovr, lcod_fle_ovr) into (name, value) [was 42340x6, now 84680x6]
#> rename: renamed one variable (overlap)
#> mutate: changed 84,680 values (100%) of 'overlap' (0 new NA)
# Since it's mainly an difference in elevation between quarters, I just do it for Q1 here
p1 <- plot_map_fc +
geom_raster(data = filter(all_pred2, quarter == 1),
aes(x = X*1000, y = Y*1000, fill = value)) +
scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
facet_grid(year ~ overlap) +
theme(legend.text = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.25, 'cm'),
plot.margin = unit(c(0, 0, 0, 0), "cm"),
plot.title = element_text(size = 10)) +
ggtitle("Quarter 1") +
NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining
p2 <- plot_map_fc +
geom_raster(data = filter(all_pred2, quarter == 4),
aes(x = X*1000, y = Y*1000, fill = value)) +
scale_fill_viridis_c(trans = "sqrt", name = "Overlap") +
facet_grid(year ~ overlap) +
theme(legend.text = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.25, 'cm'),
plot.margin = unit(c(0, 0, 0, 0), "cm"),
plot.title = element_text(size = 10)) +
ggtitle("Quarter 4") +
NULL
#> filter: removed 42,340 rows (50%), 42,340 rows remaining
# Over space
p5 <- all_pred %>%
distinct(year, quarter, .keep_all = TRUE) %>%
ggplot(aes(year, lcod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) +
labs(color = "Quarter", x = "Year", y = "Overlap index") +
guides(fill = "none") +
stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
geom_point() +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#> distinct: removed 592,704 rows (>99%), 56 rows remaining
p6 <- all_pred %>%
distinct(year, quarter, .keep_all = TRUE) %>%
ggplot(aes(year, scod_fle_ovr_tot, color = factor(quarter), fill = factor(quarter))) +
labs(color = "Quarter", x = "Year", y = "Overlap index") +
guides(fill = "none") +
stat_smooth(method = "gam", formula = y~s(x, k = 3)) +
geom_point() +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#> distinct: removed 592,704 rows (>99%), 56 rows remaining
((p1 | p2) / ((p5 | p6) + plot_layout(guides = "collect"))) + plot_layout(heights = c(1.1, 1)) + plot_annotation(tag_levels = "A") & theme(aspect.ratio = 1)
ggsave("figures/spatial_overlap.pdf", width = 20, height = 20, units = "cm")Conditional effects
# Scale with respect to data!
nd_depth <- data.frame(depth = seq(min(cpue2$depth), max(cpue2$depth), length.out = 100)) %>%
mutate(depth_ct = depth - mean(cpue2$depth),
depth_sq = depth_ct * depth_ct,
depth_sq_sc = (depth_sq - mean(cpue2$depth_sq)) / sd(cpue2$depth_sq),
depth_sc = (depth - mean(cpue2$depth)) / sd(cpue2$depth),
temp_sq_sc = 0,
temp_sc = 0,
oxy_sc = 0,
sal_sc = 0,
year = 1999L,
fyear = as.factor(1999),
fsubstrate = "mud")
#> mutate: new variable 'depth_ct' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with 100 unique values and 0% NA
#> new variable 'depth_sc' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with one unique value and 0% NA
#> new variable 'temp_sc' (double) with one unique value and 0% NA
#> new variable 'oxy_sc' (double) with one unique value and 0% NA
#> new variable 'sal_sc' (double) with one unique value and 0% NA
#> new variable 'year' (integer) with one unique value and 0% NA
#> new variable 'fyear' (factor) with one unique value and 0% NA
#> new variable 'fsubstrate' (character) with one unique value and 0% NA
# Q1
margin_small_cod_pred_q1_depth <- predict(mcod_s_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q1_depth <- predict(mcod_l_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q1_depth <- predict(mfle_q1_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
margin_small_cod_pred_q4_depth <- predict(mcod_s_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q4_depth <- predict(mcod_l_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q4_depth <- predict(mfle_q4_space,
newdata = nd_depth,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margins_depth <- bind_rows(margin_small_cod_pred_q1_depth,
margin_large_cod_pred_q1_depth,
margin_fle_pred_q1_depth,
margin_small_cod_pred_q4_depth,
margin_large_cod_pred_q4_depth,
margin_fle_pred_q4_depth) %>%
mutate(species = ifelse(species == "flounder", "Flounder", species),
species = ifelse(species == "small_cod", "Small cod", species),
species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)
# Temperature
nd_temp <- data.frame(temp = seq(min(cpue2$temp), max(cpue2$temp), length.out = 100)) %>%
mutate(temp_ct = temp - mean(cpue2$temp),
temp_sq = temp_ct * temp_ct,
temp_sq_sc = (temp_sq - mean(cpue2$temp_sq)) / sd(cpue2$temp_sq),
temp_sc = (temp - mean(cpue2$temp)) / sd(cpue2$temp),
depth_sq_sc = 0,
depth_sc = 0,
oxy_sc = 0,
sal_sc = 0,
year = 1999L,
fyear = as.factor(1999),
fsubstrate = "mud")
#> mutate: new variable 'temp_ct' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq' (double) with 100 unique values and 0% NA
#> new variable 'temp_sq_sc' (double) with 100 unique values and 0% NA
#> new variable 'temp_sc' (double) with 100 unique values and 0% NA
#> new variable 'depth_sq_sc' (double) with one unique value and 0% NA
#> new variable 'depth_sc' (double) with one unique value and 0% NA
#> new variable 'oxy_sc' (double) with one unique value and 0% NA
#> new variable 'sal_sc' (double) with one unique value and 0% NA
#> new variable 'year' (integer) with one unique value and 0% NA
#> new variable 'fyear' (factor) with one unique value and 0% NA
#> new variable 'fsubstrate' (character) with one unique value and 0% NA
# Q1
margin_small_cod_pred_q1_temp <- predict(mcod_s_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q1_temp <- predict(mcod_l_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species ="large_cod", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q1_temp <- predict(mfle_q1_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 1)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
# Q4
margin_small_cod_pred_q4_temp <- predict(mcod_s_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "small_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_large_cod_pred_q4_temp <- predict(mcod_l_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "large_cod", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margin_fle_pred_q4_temp <- predict(mfle_q4_space,
newdata = nd_temp,
se_fit = TRUE,
re_form = NA) %>%
mutate(species = "flounder", quarter = 4)
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'quarter' (double) with one unique value and 0% NA
margins_temp <- bind_rows(margin_small_cod_pred_q1_temp,
margin_large_cod_pred_q1_temp,
margin_fle_pred_q1_temp,
margin_small_cod_pred_q4_temp,
margin_large_cod_pred_q4_temp,
margin_fle_pred_q4_temp) %>%
mutate(species = ifelse(species == "flounder", "Flounder", species),
species = ifelse(species == "small_cod", "Small cod", species),
species = ifelse(species == "large_cod", "Large cod", species))
#> mutate: changed 600 values (100%) of 'species' (0 new NA)
# Plot!
p1 <- ggplot(margins_depth,
aes(depth, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
color = species, fill = species)) +
geom_line() +
facet_wrap(~quarter, scales = "free") +
geom_ribbon(alpha = 0.4, color = NA) +
coord_cartesian(xlim = c(10, 170)) +
theme(aspect.ratio = 1) +
labs(x = "Depth (m)", y = "Biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
p2 <- ggplot(margins_temp,
aes(temp, exp(est), ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se),
color = species, fill = species)) +
geom_line() +
facet_wrap(~quarter, scales = "free") +
geom_ribbon(alpha = 0.4, color = NA) +
coord_cartesian(xlim = c(1, 14)) +
theme(aspect.ratio = 1) +
labs(x = "Temperature (°C)", y = "Biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
(p1 / p2) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")
# ggsave("figures/supp/conditional_effects.pdf", width = 20, height = 20, units = "cm#")
# Scaled
p3 <- margins_temp %>%
ungroup() %>%
group_by(species, quarter) %>%
mutate(max_est = max(exp(est)),
est_sc = exp(est) / max_est) %>%
ggplot(aes(temp, est_sc, color = species, fill = species)) +
geom_line(size = 1.3) +
facet_wrap(~quarter, scales = "free") +
coord_cartesian(xlim = c(1, 14)) +
theme(aspect.ratio = 1) +
labs(x = "Temperature (°C)", y = "Scaled biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#> new variable 'est_sc' (double) with 595 unique values and 0% NA
p4 <- margins_depth %>%
ungroup() %>%
group_by(species, quarter) %>%
mutate(max_est = max(exp(est)),
est_sc = exp(est) / max_est) %>%
ggplot(aes(depth, est_sc, color = species, fill = species)) +
geom_line(size = 1.3) +
facet_wrap(~quarter, scales = "free") +
coord_cartesian(xlim = c(10, 170)) +
theme(aspect.ratio = 1) +
labs(x = "Depth (m)", y = "Scaled biomass density", color = "", fill = "") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (species, quarter)
#> mutate (grouped): new variable 'max_est' (double) with 6 unique values and 0% NA
#> new variable 'est_sc' (double) with 595 unique values and 0% NA
(p3 / p4) + plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")
# ggsave("figures/conditional_effects.pdf", width = 20, height = 20, units = "cm")Weighted quantiles
head(all_pred) %>% as.data.frame()
#> X Y depth year lon lat substrate quarter oxy temp
#> 1 450 5984 11 1993 14.23718 54.00188 sand 1 8.713132 2.821855
#> 2 454 5984 11 1993 14.29820 54.00225 sand 1 8.743857 2.811376
#> 3 458 5984 11 1993 14.35922 54.00259 sand 1 8.745697 2.806495
#> 4 462 5984 12 1993 14.42024 54.00290 sand 1 8.759304 2.802808
#> 5 466 5984 12 1993 14.48127 54.00318 sand 1 8.765212 2.805644
#> 6 470 5984 11 1993 14.54229 54.00343 sand 1 8.805174 2.792613
#> sal sub_div ices_rect density_saduria biomass_spr biomass_her
#> 1 7.293935 24 37G4 0 NA NA
#> 2 7.346682 24 37G4 0 NA NA
#> 3 7.383028 24 37G4 0 NA NA
#> 4 7.418652 24 37G4 0 NA NA
#> 5 7.501428 24 37G4 0 NA NA
#> 6 7.585449 24 37G4 0 NA NA
#> biomass_spr_sd biomass_her_sd depth_ct depth_sq depth_sq_sc depth_sc
#> 1 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 2 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 3 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> 4 NA NA -35.63035 1269.522 0.4131041 -1.306567
#> 5 NA NA -35.63035 1269.522 0.4131041 -1.306567
#> 6 NA NA -36.63035 1341.783 0.4698617 -1.343238
#> temp_ct temp_sq temp_sq_sc temp_sc oxy_sc sal_sc fyear fsubstrate
#> 1 -3.484890 12.14446 0.5884181 -1.305532 1.178372 -1.098960 1993 sand
#> 2 -3.495369 12.21760 0.5969916 -1.309458 1.190285 -1.082916 1993 sand
#> 3 -3.500250 12.25175 0.6009941 -1.311287 1.190999 -1.071861 1993 sand
#> 4 -3.503937 12.27757 0.6040213 -1.312668 1.196275 -1.061026 1993 sand
#> 5 -3.501101 12.25771 0.6016927 -1.311606 1.198566 -1.035849 1993 sand
#> 6 -3.514132 12.34912 0.6124078 -1.316487 1.214061 -1.010293 1993 sand
#> est est_non_rf est_rf omega_s epsilon_st small_cod large_cod
#> 1 -3.128299 -2.913320 -0.2149794 -0.3565650 0.1415856 0.04379222 6.078414
#> 2 -3.475978 -2.915583 -0.5603953 -0.7218136 0.1614182 0.03093157 5.816645
#> 3 -3.803139 -2.919747 -0.8833924 -1.0647534 0.1813610 0.02230065 5.483502
#> 4 -3.886795 -2.811876 -1.0749196 -1.2701738 0.1952542 0.02051097 5.240582
#> 5 -4.076864 -2.810417 -1.2664467 -1.4755941 0.2091473 0.01696057 4.760898
#> 6 -4.302368 -2.921767 -1.3806012 -1.6096988 0.2290976 0.01353646 4.322854
#> fle scod_fle_ovr scod_fle_ovr_tot lcod_fle_ovr lcod_fle_ovr_tot
#> 1 13.58323 6.025629e-08 0.3323209 9.576900e-07 0.5491266
#> 2 13.07450 4.096655e-08 0.3323209 8.821231e-07 0.5491266
#> 3 12.76604 2.883873e-08 0.3323209 8.119809e-07 0.5491266
#> 4 13.48553 2.801926e-08 0.3323209 8.197456e-07 0.5491266
#> 5 14.17812 2.435911e-08 0.3323209 7.829590e-07 0.5491266
#> 6 14.43847 1.979833e-08 0.3323209 7.239743e-07 0.5491266
all_pred_sub <- all_pred %>%
dplyr::select(small_cod, large_cod, fle, year, temp, oxy, sal, depth, quarter) %>%
pivot_longer(c(small_cod, large_cod, fle)) %>%
rename(Species = name, density = value)
#> pivot_longer: reorganized (small_cod, large_cod, fle) into (name, value) [was 592760x9, now 1778280x8]
#> rename: renamed 2 variables (Species, density)
wm <- all_pred_sub %>%
pivot_longer(c(temp, oxy, sal, depth)) %>%
mutate(name = ifelse(name == "temp", "Temperature (°C)", name),
name = ifelse(name == "sal", "Salinity", name),
name = ifelse(name == "depth", "Depth", name),
name = ifelse(name == "oxy", "Oxygen (ml/L)", name)) %>%
group_by(year, quarter, Species, name) %>%
summarise("Weighted median" = hutils::weighted_quantile(v = value, w = density, p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = value, w = density, p = c(0.75))) %>%
rename(env_var = name) %>%
mutate(Species = ifelse(Species == "fle", "Flounder", Species),
Species = ifelse(Species == "large_cod", "Large cod", Species),
Species = ifelse(Species == "small_cod", "Small cod", Species))
#> pivot_longer: reorganized (temp, oxy, sal, depth) into (name, value) [was 1778280x8, now 7113120x6]
#> mutate: changed 7,113,120 values (100%) of 'name' (0 new NA)
#> group_by: 4 grouping variables (year, quarter, Species, name)
#> summarise: now 672 rows and 7 columns, 3 group variables remaining (year, quarter, Species)
#> rename: renamed one variable (env_var)
#> mutate (grouped): changed 672 values (100%) of 'Species' (0 new NA)
ggplot(wm, aes(year, `Weighted median`, color = Species, fill = Species)) +
geom_line() +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
ggplot(wm, aes(year, `Weighted median`, color = Species)) +
geom_line(size = 0.8) +
labs(linetype = "Quarter", x = "Year") +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#ggsave("figures/weighted_quantiles.pdf", width = 20, height = 20, units = "cm")
# Scaled
wm %>%
group_by(Species, quarter, env_var) %>%
mutate(`Weighted median scaled` = `Weighted median` - mean(`Weighted median`)) %>%
ggplot(aes(year, `Weighted median scaled`, color = Species)) +
geom_line(size = 0.8) +
labs(linetype = "Quarter", x = "Year") +
ggh4x::facet_grid2(env_var ~ quarter, scales = "free_y", independent = "y") +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#> group_by: 3 grouping variables (Species, quarter, env_var)
#> mutate (grouped): new variable 'Weighted median scaled' (double) with 656 unique values and 0% NA
#ggsave("figures/supp/weighted_quantiles_scaled.pdf", width = 20, height = 20, units = "cm")
# ggplot(wm, aes(year, Median, color = Species, fill = Species)) +
# geom_ribbon(aes(ymin = `1st quartile`, ymax = `3rd quartile`), alpha = 0.4, color = NA) +
# facet_wrap(env_var ~ quarter, scales = "free", ncol = 2) +
# scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
# scale_fill_brewer(palette = "Paired", name = "Species", direction = -1)
#
# ggsave("figures/supp/weighted_quantiles_range.pdf", width = 20, height = 20, units = "cm")[Extra for bentfish - calculate biomass index per subdivision]
Compare with Q4 from previous index
cpue3 <- cpue2 %>%
mutate(tot_cod = large_cod + small_cod) %>%
filter(quarter == 4)
#> mutate: new variable 'tot_cod' (double) with 7,870 unique values and 0% NA
#> filter: removed 5,271 rows (59%), 3,595 rows remaining
spde <- make_mesh(cpue3, xy_cols = c("X", "Y"),
n_knots = 200,
type = "kmeans", seed = 42)
cod_ind <- sdmTMB(tot_cod ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = cpue3,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
fle_ind <- sdmTMB(flounder ~ 0 + fsubstrate + fyear + depth_sc + depth_sq_sc +
temp_sc + temp_sq_sc + oxy_sc + sal_sc,
data = cpue3,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1))
sanity(cod_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
sanity(fle_ind)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large# cod_ind2 <- sdmTMB(tot_cod ~ 0 + fyear + s(depth_sc),
# data = cpue3,
# mesh = spde,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# fle_ind2 <- sdmTMB(flounder ~ 0 + fyear + s(depth_sc),
# data = cpue3,
# mesh = spde,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# sanity(cod_ind2)
# sanity(fle_ind2)Now try old data as well
# d_old <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
#
# # Filter to match the data I want to predict on
# d_old <- d_old %>%
# filter(year > 1992 & year < 2020 & quarter == 4)
#
# # Calculate standardized variables
# d_old <- d_old %>%
# mutate(depth_sc = (depth - mean(depth))/sd(depth),
# X = X/1000,
# Y = Y/1000,
# year = as.integer(year),
# fyear = as.factor(year)) %>%
# drop_na(depth) %>%
# rename("density_cod" = "density") # to fit better with how flounder is named
#
# nrow(d_old)
# nrow(cpue3)
#
# # Plot comparison
# p1 <- ggplot() +
# geom_histogram(data = d_old, aes(density_cod, fill = "old"), alpha = 0.5) +
# geom_histogram(data = cpue3, aes(tot_cod, fill = "new"), alpha = 0.5) +
# scale_y_continuous(trans = "log")
#
# p2 <- ggplot() +
# geom_histogram(data = d_old, aes(density_fle, fill = "old"), alpha = 0.5) +
# geom_histogram(data = cpue3, aes(flounder, fill = "new"), alpha = 0.5) +
# scale_y_continuous(trans = "log")
#
# p1 + p2# spde_old <- make_mesh(d_old, xy_cols = c("X", "Y"),
# n_knots = 200,
# type = "kmeans", seed = 42)
#
# cod_ind3 <- sdmTMB(density_cod ~ 0 + fyear + s(depth_sc),
# data = d_old,
# mesh = spde_old,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# fle_ind3 <- sdmTMB(density_fle ~ 0 + fyear + s(depth_sc),
# data = d_old,
# mesh = spde_old,
# family = tweedie(link = "log"),
# spatiotemporal = "AR1",
# spatial = "on",
# time = "year",
# reml = FALSE,
# control = sdmTMBcontrol(newton_loops = 1))
#
# sanity(cod_ind3)
# sanity(fle_ind3)### Cod
# SD 24
preds_cod24_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining
# SD 25
preds_cod25_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining
# SD 26
preds_cod26_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining
# SD 27
preds_cod27_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining
# SD 28
preds_cod28_sim <- predict(cod_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining
# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(4*4, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(4*4, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(4*4, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(4*4, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(4*4, nrow(preds_cod28_sim)))
# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
### Flounder
# SD 24
preds_fle24_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 513,744 rows (87%), 79,016 rows remaining
# SD 25
preds_fle25_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 446,040 rows (75%), 146,720 rows remaining
# SD 26
preds_fle26_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 457,856 rows (77%), 134,904 rows remaining
# SD 27
preds_fle27_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 498,120 rows (84%), 94,640 rows remaining
# SD 28
preds_fle28_sim <- predict(fle_ind, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 455,280 rows (77%), 137,480 rows remaining
# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(4*4, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(4*4, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(4*4, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(4*4, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(4*4, nrow(preds_fle28_sim)))
# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
div_index_sim <- bind_rows(index24_cod_sim,
index25_cod_sim,
index26_cod_sim,
index27_cod_sim,
index28_cod_sim,
index24_fle_sim,
index25_fle_sim,
index26_fle_sim,
index27_fle_sim,
index28_fle_sim) %>%
mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes
#> mutate: new variable 'est_t' (double) with 280 unique values and 0% NA
#> new variable 'lwr_t' (double) with 280 unique values and 0% NA
#> new variable 'upr_t' (double) with 280 unique values and 0% NAPlot the index and save as csv
# # First compare with previous index:
# cod_fle_index_old <- read.csv("/Users/maxlindmark/Dropbox/Max work/R/cod_condition/output/cod_fle_index.csv") %>%
# mutate(sub_div = as.character(sub_div)) %>%
# filter(sub_div %in% c("25", "28"))
#
# index_comp <- bind_rows(cod_fle_index_old %>% mutate(source = "cod_condition"),
# div_index_sim %>% mutate(source = "cod_interactions"))
#
# ggplot(index_comp, aes(year, est_t, color = source, fill = source)) +
# facet_wrap(sub_div ~ species, scales = "free") +
# geom_line() +
# geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t), alpha = 0.2, color = NA)
#
#
# ## If it isn't the same with the data, plot pred grids against each other
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2, color = NA) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
div_index_sim %>%
filter(sub_div %in% c(25, 28)) %>%
ggplot(aes(year, est_t/1000, color = species, fill = species)) +
geom_line() +
facet_wrap(~sub_div, scales = "free", ncol = 1) +
geom_ribbon(aes(year, ymin = lwr_t/1000, ymax = upr_t/1000), alpha = 0.2, color = NA) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "Species") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
#> filter: removed 168 rows (60%), 112 rows remaining
ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
geom_line() +
facet_wrap(~species, scales = "free") +
#geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
ggplot(div_index_sim, aes(year, est_t, color = species)) +
geom_line() +
facet_wrap(~sub_div, scales = "free") +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
write.csv(div_index_sim, "output/cod_fle_index.csv")